The power of random measurements: measuring Trp" on single copies of p 
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While it is known that Trp™ can be measured directly (i.e., without first reconstructing the 
density matrix) by performing joint measurements on n copies of the same state p, it is shown here 
that random measurements on single copies suffice, too. Averaging over the random measurements 
directly yields estimates of Trp™, even when it is not known what measurements were actually 
performed (so that p cannot be reconstructed). 



The standard textbook quantum measurement of an 
observable O on a given quantum system produces an 
estimate of the expectation value Tr(pO), where p is the 
density matrix of the system. This expectation value is 
linear in p. As is well-known by now [IH5] , nonlinear func- 
tions of the density matrix p, such as the purity p2 = Trp 2 
and its cousins p„ = Trp n for n > 2, can be measured di- 
rectly, too, without first having to reconstruct the whole 
density matrix. For this direct measurement method to 
work one needs n quantum systems that are all in the 
same state p, plus the ability to perform the appropriate 
joint measurement (s) on those multiple copies. 

Here we point out that estimates of the same nonlinear 
quantities can be obtained from random measurements 
on single copies as well. A random measurement can 
be assumed to be implemented by performing a random 
unitary rotation on the single copy (possibly including 
an ancilla which starts off in a standard state), followed 
by a fixed measurement on the single copy (and possibly 
on the ancilla). By averaging the measurement results 
over the random unitaries, one can directly infer esti- 
mates of Trp" [with n = 2...M, with M the Hilbert 
space dimension of the system of interest] , without hav- 
ing to reconstruct the density matrix. One point of the 
averaging procedure is that one does not have to know 
which random unitaries were in fact applied, and as a 
consequence one cannot reconstruct the density matrix 
in that case. An example of a random measurement is 
furnished by intensity measurements of speckle patterns 
resulting from light (be it two photons, or a single pho- 
ton, or a coherent laser beam) propagating through a 
disordered medium [6l [7] , and in that case the purity P2 
can (and was indeed) inferred directly from those mea- 
surements (see also [T?]). 

There is an important difference between the known 
direct method and the current random method in what 
quantity exactly is estimated. Suppose one's source does 
not produce the same state every single time, but instead 
a state pj at try j. In this case standard quantum mea- 
surements of a given observable on J instances j = 1 . . . J 
can still be described by a single density matrix, namely, 
the mean p = Y] ^ pj/ J. Since the random method only 



involves measurements on single copies, it produces, like- 
wise, an estimate of Tr(p n ). This requires no assumption 
about the quantum systems being uncorrelated or unen- 
tangled with each other, since pj is obtained by tracing 
out all degrees of freedom except those of system j. 

On the other hand, a direct measurement would 
yield an estimate of Ti(Spjj+i...j+ n -i) instead, where 
Pj,j+i—j+n-i is the joint density matrix of n systems 
j,j + l,...,j+n— 1, and S is the cyclical shift operator, 
which acts on the basis states of the n quantum systems 

as Sl^l^+l) ■ • ■ |^j+n-l) = \ipj+i)\ipj+2) ■ ■ ■ Wj)- It is 
only under the assumption that the states of the n sys- 
tems are identical and independent (i.i.) that the direct 
measurement yields Trp™. In fact, the direct measure- 
ment is eminently suited for detecting that the states are 
not identical [9 . Although the assumption of i. i. states is 
standard, it is only recently that precise conditions have 
been stated under which the approximate i.i. character 
can be inferred |10j . The required permutation invari- 
ance is easily enforced when performing measurements 
on single copies, but not when performing joint measure- 
ments on multiple copies Avoiding this difficulty is 
the main advantage of the random method. 

An N x N random unitary matrix, distributed accord- 
ing to the Haar measure, can be easily constructed by the 
method presented in |12j . One first constructs a matrix 
whose elements are independent complex Gaussian vari- 
ables, and one then performs an orthogonalization of the 
resulting random matrix (where one small pitfall needs to 
be avoided |12j). We first consider approximate results 
for random unitaries, because the resulting expressions 
are quite simple, and subsequently we will give the more 
involved exact results. 

If we consider an arbitrary submatrix V (of size M) of 
U (of size N) , with M <C N Q2] , then the real and imag- 
inary parts of its matrix elements can still be very well 
approximated by independent and normally distributed 
numbers if N is large. With this Gaussian approximation 
we can compute the following averages (we indicate av- 
erages over the distribution of random unitaries by (.)): 
first, we have 

(VkiKJ = S km S ln /N. (1) 
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Here and in all of the following we assume we have picked 
some basis {]&)}, and we write all matrix elements w.r.t. 
that basis. The normalization factor 1/N follows imme- 
diately from the fact that U, of which V is a submatrix, 
is unitary, so that Yli=i UkiU* nl = 5km- Higher-order 
averages follow from the Isserlis ("Gaussian-moments") 
theorem [14j . In particular, the only nonzero averages 
arise from products of 2K factors of the form 
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We now apply the preceding approximate results to the 
following scenario. Consider an "input" density matrix 
p ln of size M x M. Embed the system in a larger Hilbcrt 
space of size N, by constructing a new N x N density 
matrix by adding zero matrix elements. Then apply a 
random unitary U to the larger matrix. Finally, consider 
measurements in a fixed M-dimensional (sub) basis 
The probability Prob(fc) of finding measurement outcome 
k is given by 



Prob(/c) 
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This expectation value depends on what V is, of course, 
but its average is given simply by 



(Prob(fc)) = J^P^n (V mk V: k ) = 1/N, 



(4) 



where we used that Tr(p ln ) = J2 m Pn 



1 . Defining 



Pn{k) = (Prob(fc)™), the following averages are obtained 
by using the Isserlis theorem (up to order n = 4; sub- 
sequent orders can be easily obtained, too, but for our 
purposes this will do) 



P2(k) 
Pi{k) 
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P2]/N 2 , 
-3p 2 + 2p 3 }/N 3 , 



[l + 3p 2 2 + 6 P2 + S P3 + 6p 4 ] /N 4 
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(5b) 
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where we defined p n — Tr ((/9 ln )™)- Inverting these equa- 
tions gives estimates of p n in terms of the measurable 
quantities on the left-hand sides. We denote those esti- 
mates by an overbar, e.g., p2 — N 2 P2(k) — 1. We refrain 
from giving the other inverse relations now, as we will 
give the exact relations below in ([9]). 

We can also compute standard deviations in the 
(mean) estimates. For example, assuming we average 
the results for one value of k over 2V ran d random uni- 
taries, then the statistical error in the estimate of the 
purity is 

V^rand - lA(p 2 ) = ^ 4 P 2 + 2p 2 + 8p 3 + 6 Pi . (6) 

This is an increasing function of P2,Pz,Pi, so that the 
variance is largest for a pure state and smallest for the 
totally mixed state p in = t/M. 



In an actual experiment one may not know exactly 
what the values of N and/or M are (for instance, this 
is the case in the speckle experiments of Refs. [51 [7]). 
In such a case N can be directly estimated from P\{k) 
through N — 1/P\(k). So, we would use 



P2 = 



pm 
aw 2 



- 1, 



(7) 



instead (such estimates we indicate by a tilde) . Now this 
estimate p2 has a smaller variance than p 2 has, simply 
because the errors in P\{k) and P2{k) are positively cor- 
related. It is, therefore, better to use p~2 as estimate for 
P2, even when N is in principle known. The numerical 
results given below will confirm this, also for the exact 
result for p2- For the estimates ps and p4, however, there 
is not much difference between the two methods. 

When N is not very large, equations ^ and hence ^ 
are not correct. The exact results, which can be extracted 
from Refs. [TS] and [TB], are still given by ^ upon mul- 
tiplication of P n (k) by the correction factor C„, where 



(1 + 1/N)(1 + 2/N) . . . (1 + (n - 1)/N). (8) 



Note that these factors depend only on N, not on M, 
and the results are valid even when M — N. This then 
leads to the inverse formulas: 



P2 
P3 
Pi 



D 2 P 2 (k) - 1, 
l[D 3 P 3 (k)-l-3p 2 ], 



l[D 4 P 4 (k)-l 
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with D n = (N + n - 1)1/ (N - 1)1 Taking into account 
the correction factors Q leads to different values for the 
statistical errors in estimates. It is still true that pure 
states lead to the largest errors; for those we get 



VXand - lA(pa) = 



' 24(l + l/iV) 
(l + 2/JV)(l + 3/JV) 



-4. (10) 



The right-hand side (slowly) increases with increasing N, 
from v/52/7 for N = 4 to V^O for N -> oo. 

In order to illustrate the method and the meanings 
of N and M, we consider the following examples here: 
(i) Suppose we have a single photon occupying one of 
M input modes. We then apply a random linear optics 
transformation that involves N — M ancilla modes. The 
photon now ends up being coherently distributed over N 
output modes. We then estimate the probability Prob(/c) 
with which the photon ends up in one of a fixed set of M 
output modes fc = 1 . . . M. This is an example akin to 
that considered in [SJ [7] ■ 

(ii) Suppose our system of interest consists of 2 qubits, 
so that M = 4. Suppose we have an ancilla qubit in a 
fixed state |0), and we apply a random unitary operation 
to the 3 qubits. In this case, N — 8. We then perform 
measurements on each of the three qubits separately in 
the standard basis. We measure the probability Prob(/c) 
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of the two qubits ending up in one of the M = 4 combi- 
nations k — 00, 01, 10, 11 and the ancilla ending up in |0) 
(thus measuring only a M-dimensional subspace). 

(ii') There is no need for any ancillas if dealing with a 
fixed and known number of qubits, say Q. In that case, 
we simply have N = M = 2® . We consider only case 
(ii') in the following numerical results. 

We assume that we run an experiment with a fixed 
random ( "unknown" ) unitary of size N sufficiently many 
times that we get a very good estimate of Prob(fc) for 
each k for the given unitary and the given input state (of 
size M). Subsequently we average over iV ran d random 
unitaries to obtain P n (k) = (Prob(/c) n ). From those re- 
sults we estimate the values of P2,P3,P4- The first ex- 
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FIG. 2: Same as previous Figure, but using the estimate N as 
1/Pi{k) in (£}. Here we have Ap 2 = 0.09, Ap 3 = 0.17, Ap 4 = 
0.27. The mean estimates obtained from pooling all data 
(which are the same "raw" data as in Fig. [TJ from the 100 
trials for p„ are p 2 = 1.005 [which is indeed better than pa], 
p~3 = 1.01 and p~4 — 1.02, all consistent with the statistical 
errors in the mean (which are 10 times smaller than Ap n ). 



independent (and known) measurements. Here we show 
that with just a moderate overhead one can obtain good 
estimates of p2,P3,P4 for generic (i.e. randomly picked 
[15] ) states. In Fig. [3] results are displayed for 200 generic 



FIG. 1: This plot shows, for a pure two-qubit state, the es- 
timated values of p2,P3,P4 (blue: P2, red: p3, green: p±) for 
100 trials, each trial using just one value of k, containing an 
average over -/V ran d = 100 random unitaries, and using N = 4 
m |9f. The mean standard deviations (over 100 trials) were 
Ap2 = 0.282 [note that this agrees with the result (10 1, since 
^52/7/99 « 0.274], Ap 3 = 0.21, Ap 2 = 0.29. The mean esti- 
mates obtained by pooling all data from the 100 trials for p n 
are: p 2 = 0.990, p 3 = 1.01 and p4 = 1.02, which are all con- 
sistent with their mean standard deviations (10 times smaller 
than the Ap n given above). 

ample we consider corresponds to case (ii') mentioned 
above, where we have two qubits. In Figs. [T] and [2] we 
plot results for pure input states, where we use the re- 
sults for just 1 value of k to estimate p n , in two different 
ways: using the exact value N = 4 (Fig. [T]) or using the 
estimate N « 1/ (Prob(fc)) (Fig. [jjj). The results show 
how the latter method is more accurate for estimating 
purity. The same data are used in the two Figures, so 
that all differences between them are entirely due to the 
different analysis of those data. This different analysis 
reduces the statistical variation in p2 , but not in and 
p4.. In addition, the plots show that the statistical errors 
in p 2 ,P3,P4 are strongly correlated in the latter case. 

In the remaining figures we perform an additional av- 
erage over the M different values of k, leading to smaller 
(by a factor of about \fM) error bars. 

Performing tomography on two qubits would require 15 




FIG. 3: Scatter plot of estimated values of P2 (blue crosses), 
P3 (red circles), and p4 (green diamonds), versus their actual 
values for 200 randomly picked two-qubit input states. Here 
an average is taken over iV ran( i = 30 random unitaries, as 
well as over 4 measurement outcomes. For convenience, the 
dashed line gives the diagonal on which estimated and actual 
values agree. 

two-qubit states, using AT ran d = 30. 

In Fig. [2] we show (for five qubits) that the number of 
random unitaries needed to obtain a fixed-size error bar 
does not increase with the number of qubits. For iV ran d = 
30 one still obtains good estimates: in fact, the error bars 
decrease (roughly as l/yM) when going to more and 
more qubits, just because the number M of measurement 
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results one can average over increases exponentially with 
the number of qubits, while the variance ( 10 1 increases 
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FIG. 4: Same as Fig. [3] but for randomly picked five-qabit 
states (M = N = 32). Averaging over the same number 
of random unitaries (here iV ran d = 30) produces a smaller 
statistical error for larger systems. 
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FIG. 5: The standard deviation in the mean estimates for 
f>2 (blue crosses), p 3 (red circles), and P4 (green diamonds) 
after averaging over N la ,nd = 30 random unitaries, for pure 
multi-qubit states, as a function of the number of qubits. 



only very slowly. This is illustrated for pure multi-qubit 
states in Fig. [5j It shows that the statistical error in 
the estimate of Tip" for n = 2,3,4 first increases with 
the number of qubits before (at > n qubits) it starts to 
decrease monotonically. 

In conclusion then, using the ideas of random matrix 
theory, we showed that nonlinear functions of the density 
matrix such as Tvp n can be directly obtained from appro- 
priately averaged random measurements on single copies. 
No assumptions are needed on the independence of the 
copies, nor on their states being identical. This contrasts 
the random method with so-called direct measurements 
on n identical copies QJ[S]. 

Moreover, one does not need to know which random 
measurements were actually performed, because the av- 
eraging procedure keeps all information about the eigen- 
values of p, which is all that is needed to estimate Trp n . 
One does need to verify that the random unitaries have 
been drawn from the appropriate ensemble. There are 
two tests one could perform: first of all, the definition of 
the ensemble is that it is unitarily invariant. This means 
in our context that all averages (Prob(/c) n ) should be 
independent of k. This is a statistically testable prop- 
erty. In addition, one can apply the random measure- 
ments to known input states, so that the values of those 
fc-independent averages are known. 

Importantly, the number of unitaries over which one 
has to average in order to obtain a fixed error bar in the 
estimates of Trp n scales very favorably with the Hilbert 
space dimension of one's system: in fact, this number 
even tends to decrease. For two qubits this amounts to 
needing a small overhead as compared to full quantum- 
state tomography, but for larger systems (more than, say, 
four qubits) the random method requires (far) fewer re- 
sources than does full quantum-state tomography. 
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